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I. INTRODUCTION 

The design and selection of a nozzle for a hypersonic scramjet must be based 
on a compromise between internal and external flow requirements related to 
vehicle lift, drag, pitching moment, thrust, structural and weight limita- 
tions. The design process involves a complex study based on engineering 
analysis and refinements using complex computer programs. A logical design 
sequence consists of first obtaining a satisfactory range of aerodynamic para- 
meters utilizing simplified analysis * ’ and then narrowing the range of para- 
meters through more accurate but complex calculations. 

This report describes a two dimensional second-order characteristic procedure 
capable of analyzing the aerodynamic performance of typical nozzle configura- 
tions selected from simplified analysis as shown in Figure (1). 



FIGURE 1. - TYPICAL SCRAMJET NOZZLE 

However, the calculation procedure is not limited to these configurations but 
can be readily adapted to calculate other two dimensional configurations. This 
generality results from the use of three coordinates systems, axisymmetric, 
axially expanding (source type flow) and Cartesian (plane two dimensional). 
Automatic provisions for switching from axially expanding to Cartesian coordi- 
nates at a specified axial station and multiple source origins are provided 
for as a user option. This unique feature allows the lateral nozzle area varia- 
tion, as in Figure (2), to be accounted for in a quasi-two dimensional manner. 



TR 186 


Page 2 


A higher order calculation would involve a fully three dimensional calcu 
lation which would locate the lateral waves. 



FIGURE 2. MULTIPLE SOURCE ORIGINS FOR LATERAL 
NOZZLE AREA VARIATION. 

The working fluid is assumed to be a hydrogen-air mixture in frozen or 
chemical equilibrium. The mixture thermodynamics is expressed via curve 
fits; i.e., individual species curve fits for frozen flow^ and mixture fits 

5 

for equilibrium flow . 

The following boundary conditions are provided for in the calculation. 

(1) Wall boundaries 

(2) Shock boundaries (Equilibrium flow only) 

(3) Contact surface (Equilibrium with shock only) 

(4) Underexpansion interaction (Equilibrium only) 

(5) Overexpansion interaction (Equilibrium only) 

(6) Prandtl-Meyer (Equilibrium only) 

While the nozzle may be over or underexpanded at the cowl, as a user option 

no external interaction need be selected. It is always assumed that the 

nozzle has a centerbody or lower wall, 
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As part of the calculation a running integration of pressure on the nozzle 
surfaces is performed which yields thrust, lift and pitching moment in 
vehicle coordinates. The pitching moment calculation requires the location 
of the moment axis to be specified. Appendix III illustrates the method. 

Section II (a,b) describes the basic flow equations for a rotational non- 
homoentropic gas-mixture. The derivation of the characteristic equations is 
given in Appendix I and the thermodynamic curve fit data is given in Appendix 
II. Section 11(c) describes the numerical scheme and grid employed, while 
Section III discusses the various boundary conditions. Some sample calcula- 
tions are presented in Section IV, and Section VI contains concluding remarks. 
Reference (7) contains a description of the program and a sample input. 
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H. BASIC EQUATIONS 

A. Frozen Chemistry - The equations governing the two dimensional, 
axi symmetric, or axially expanding inviscid flow of a gas mixture, with 
frozen chemistry, may be written as follows: 

Continuity: 

S-Momentum: 

N-Momentum: 

Energy: 

Species Con- 
servation: 

State: 

where 0^ and Jg are the axisymmetric and axially expanding (source) terms re- 
spectively. By straightforward algebraic maniupulation , the above equations 
may be cast into characteristic form (as done in Appendix I). 


L CpJI . ) . ± pq liL + j £3. sine + ^ cose - 0 

an i y 


as 


pq 


2 x 


la. + ir = o 

as as 


pq 2 86 + = 0 

9s an 


pq 


(Y -1) 


2 “p ^ as 


21- n i£ = 


9S 


( 1 ) 

( 2 ) 

(3) 

(4a) 


co co 


3a . 

^1=0 (1=1, NSP) 


n _ £X 

p w 


y M‘ 

1 rv\ r 


(5a) 

(6a) 


Let C, denote an up-running and C_ denote a down-running characteristic. Then, 
along a C + characteristic, whose slope is expressed by: 

3 x = tan (e±p f ) (7a) 

the compatibility relation may be written 


Siny^ COSy. 



d In p ± de+ (J^ 


sine i_ , coss\ 
y u 2 x ' 


siny f 

cos(e±y f ) 


dx = 0* 


(8a) 


* 


The use of d(LnP) in place of — jj considerably improves the accuracy of the 
results for a given mesh spacing. 
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It is to be noted that at a point in the flow (x,y), the properties are 
completely specified by q, T, p, 6 and a.. (i = l» NSP). Other variables 
may be calculated as follows. 

The molecular weight is expressed by 

w= 1 % (9) 

m.j 

hence the mixture's gas constant is 

R 3 

W. 

The density is obtained from the equation of state 


y m; 


p - 


WP _ 
T W 


The thermodynamic properties C (T), h. (T) and s^(T) are tabulated poly- 

K 'j i . 

nomials, a description of which may be found in Appendix II. 


The specific heat of the mixture is expressed by 
NSP 

C n = l C n a. 

Pf i=i p i, 1 


(10) 


and the ratio of specific heats by 

. V 

' , '~V 7 V 

The local frozen Mach number is 

k 


M f “ 


M q y R 

00 ^ * 00 OQ 


/r yR 

and the Mach angle is given by 


(11) 


(12) 


**f = si* 1 m 


(13a) 
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B. Equilibrium Chemistry - The equations governing the two dimen- 
sional, axisymmetric, or axially expanding inviscid flow of a gas mixture 
in chemical equilibrium may be written as follows: 


Continuity: 

3 / \ ^ 30 . ^ sine , ! _ cose _ n 

is (pq) + pq aF + °1 pq y + J 2 pq x ' 0 

(l) 

S-Momentum: 

nD M + l£ = o 
pq 3S as u 

(2) 

N-Momentum: 

2 39 , 3£ „ n 

pq i? + a£ ' 0 

(3) 

Conservation of ~u. 

_ .. , = 0 where H = h + i q^ 

Stagnation Enthalpy: 3s 2 

(4b) 


Constancy of Equlva 
lence Ratio Along 
Stream! ines : 


o 

9S 


(5b) 


Caloric Equation 
of State: 


p/p r « constant 


(6b) 


where the equilibrium isentropic exponent is given by 


r = f (h, p,$)* 

Then along the C + characteristic whose slope is given by 


(14) 


dx = tan(e± * , e ) 


(7b) 


the compatibility relation may be written as 


sinu cosy . smy 

— V— d in P ± de + (Jl f + 0 2 s^.) dx=0(8b) 


Thus, at a point (x,y) in the flow the properties q, h, p, p, r, 4>,e are known. 

*Fitsfor r have been e'x'pressed" "in polynomial form from properties tabulated" 
in Reference ( 5 ) as described in Appendix II. 
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The local equilibrium Mach number and Mach angle are given by 
M e = q/a e 
and u e = sin" 1 

. e 

2 

where a 0 = r P^ p 

C. Numerical Procedure and Characteristic Network - Figure 
(3) depicts the global grid ordering scheme used in the present program. 
While a free running characteristic network is used the program orders and 
stores data along down-running (C_) characteristics . The only exception 
being the initial data line which must be a non-characteristic line. The 
marching proceeds from one C_ line to another until the desired flow field 
is calculated. 

A typical characteristic mesh is depicted in Figure (4), properties being 
known along the line AB and to be determined at the point C. 



i 


(15) 
(13b) 

(16) 
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Let 

* a tan(e+y)^ + 3 tan(e+p)^ 
Mg - a tan(e-y)g + g tan(e-y)^ 
and Mg = a tan + g tan 


(17a) 

(17b) 

(17c) 


The Mach angles (y) are the local values corresponding to either frozen or 
equilibrium flow. The a and 6 in the above and following equations are 
used as artifices in averaging properties along characteristics. In a first 
approximation a would be set equal to one and g equal to zero, thus fixing 
values at the points A, B or 0. Once properties at point C are determined, 
the coefficients involved in the calculation are averaged by setting both 
a and e equal to one-half. This corresponds to the second iteration. 

Writing Equation (7) in finite difference form 


and 


y C’ y A 

— £. = M 

*C-*A 1 


y C~ y B 

— E. = M 

x c -x B ; 


Solving the above yields 


X C = 


. y B-V M l x A- M 2 x B 


M r M 2 


(18a) 


(18b) 


(19a) 


y C = y A + H 1<V X A J 


(19b) 


For frozen flow let 


A, = a( S1 ' np cosw ) + B( SlnM C0SV ) 

1 * A Y C 


(20a) 


g - /S iny COSy \ + w Siny COSy v 

1 Y Y , 


(20b) 
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for equilibrium flow the Mach number and angles y are the local equilibrium 
values and the ratio of specific heats y is replaced by the equilibrium ex- 
ponent r. Similar remarks apply to all further characteristic coefficients. 


a - i r ; 5in9 siny x . „/ Sine siny x . 
2 1 a 'y'7ToY(0+"y) a ^7 cos'(6+y)y 


(21a) 


, , r ,cose siny \ , „/Cose siny \ , 

+ J 2 [a( x cos(e+y) A + B( x cos(e+y) ^ 


(21b) 


„ _ , , /sine siny i t „/Sine siny \ , 
h - J 1 la{ y-cosTe- v )\ + B( y cosT^- 7 )^ 1 

, , r ,cose siny n ^ „/Cose siny \ ■, 

+ °2 [a( x~Yosi?-y) ~ ) B + g ^' x VoV( e -y ~ ) ~ c 

Hence, along AC (C + characteristic) 

Aj (In p c ~ In p a ) + e c -e A + Ag (xq—x^) - 0 


(21c) 

(21d) 


(22a) 


and along BC (C_ characteristic) 


B 1 P B^ “ e C +e B + B 2 ^ X C’ X B^ = 0 

Solving the above equations for In p^ yields 

In p c = (A : In p A + Bj In p B + e A -e B - (A 2 +B 2 ) x c 

Ag X A + B2 x g ) / ( A j + B^) 

hence p^ - exp (In p^) and the flow inclination is 

6 C = 0 A " A i ^ ln p C - ln P A^ “ A 2 ^ x C“ x A^ 

If either the entropy or stagnation enthalpy is not constant throughout the 
flow field, the streamlines are the third family of characteristics, whose 


(22b) 


(23) 


( 24 ) 
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slope is given by 

= tan 9 , (25) 

Then,' in difference form, 


y_c i 

x (f x 



3 


(26) 


Referring to Figure (4), the point 0 lies between points A and B and can be 
located by an iterative procedure using Equation (26). Properties at point 
D are then obtained by linearly interpolating between A and B. 


The velocity at point C is found using the S-Momentum equation in the form; 


2(p c - p D ) 

q C q D [a(pq) D +S(pq) c ] 


(27) 


If the chemistry is frozen the temperature may be obtained using the Energy 
Equation (4a) along CD. That is 



T . ( Pc- p D } 

D a ( p C _ ) +3(pC) 

P D P C 


The Species Conservation Equation (4a) yields 


( 28 ) 


a, = a. (i = 1, NSP) (29) 

X “Id 

and the remaining variables are found usingEquations (6a) through (13a). 


If the chemistry is in equilibrium the Energy Equation (4b) yields the static 
enthalpy 


h C * h D + ~ q C^ 2 

and the constancy of equivalence ratio along streamlines yields 


( 30 ) 
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the curve fits for the isentropic exponent r yield 

r c = r{p c . h c , <t c ) (14) 

and the State Equation (6b) yields the density. 

The remaining variables are obtained from Equations (13b), (15) and (16). 

The calculation is then repeated with the coefficients averaged, by setting 
a and 3 equal to 1/2. If properties change significantly between these two 
sets of calculations, this generally implies too large a mesh in this region. , 

, ! • K in i ‘ , , ; " < 
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III. BOUNDARY CALCULATIONS 

A. Upper Or Lower Wall - The nozzle wall shapes, either upper 
wall (cowl) or lower wall (vehicle undersurface) are specified by polynomials 
of the form 

yj = A j ^ x “ x ^ 2 + B j( x-x i) + c j ^ 32a ) 

where x. is the origin of the wall. A maximum of 3 wall segments are permitted, 
i.e., j « 3. The wall slope e w is given by 

tan e i( * 2A-(x-x..) + B. (32b) 

W J l J 

In Figure (5), DC is the specified upper nozzle wall. The point A lies on 



c- 


FIGURE 5. WALL POINT 

the C characteristic DA on which all properties are known and point C is 
located for a = 1 and 3 = 0 by a simultaneous solution of Equation (32a) and 


X C -X A 


a tan(e+y)^ + & tante+y)^ - 


( 33 ) 


Note that this solution involves a minor iteration since (x^, y^) and hence 
6 C = 6 w^ x C^ are not known a Priori. 
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Having located point C, is known and the compatibility relation along AC, 
i.e,, Equation (22a) determines the pressure 

e fl - e r - A 9 (x r -x«\ 

In p c » In p A + — ■ Ap c '= exp(ln p c ) (34) 

the streamline equations applied along DC for either frozen or equilibrium 
flow then determine the remaining flow variables. The process is then re- 
peated for a * p 0 * Similar remarks apply for a lower wall calcula- 
tion except the characteristic relation (Equation 22b) is applied along 
a down characteristic (C_). 


B. Shock Phenomena - The program developed has the capability 
of computing the shock strength associated with an inviscid supersonic over 
or under-expansion process, and a shock propagating into a nonuniform media. 


1. Hugoniot Relations - Assume a coordinate system 
oriented, along (t direction) and normal to the shock surface (n direction) 
as shown in Figure (6). The angle sigma a is the direction cosine of the shock 
with respect to the Cartesian direction x, and u and v are the velocity com- 

A A 

ponents in the n and t directions. 


n = - sinai + coSai 

a Jr 


t = cosai + sinai 
x y 


V = - u n + v. t = ui + vi 
y 1 * y 



(35) 

(36) 

(37) 



The Rankine Hugoniot relations for a mixture in chemical equilibrium take 
the form 


Continuity: 


PjUj = P 2 u 2 


Normal 

Momentum: 


^2 

1 U 1 


^2 

p 2 + p 2 u 2 


Tangential 

Momentum: 


Energy: 


r\j 



H = h + j = constant 


State: 


p = p(p, h, *i) 


where 


4> = constant 


Employing the jump relations for a given shock angle and upstream conditions 
requires an iteration process since the mixture is calorically imperfect. 


Let 1 designate upstream conditions and 2 downstream conditions. To solve 
the jump relations knowing conditions at 1, a value for u 2 is assumed. The 
density p 2 is computed using Equation (38), p 2 is computed using Equation 
(39) and Equation (41) yields a value for h 2< The State Equation (42) then 
yields an alternate value for the density. If this value for density does 
not agree with that calculated from continuity to within a specified toler- 
ance, a new value of u 2 is assumed and this process is repeated until con- 
vergence is achieved. 

2. Shock Point Calculation - Referring to Figure (7), 
a typical shock wave calculation is performed as follows. A value of the 
shock angle is assumed, and a simultaneous solution of the equations 
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FIGURE 7. SHOCK POINT CALCULATION 


and 


y C“ y A 


X C" X A 


y C~ y B 
V- ■ = M 0 

X--X~ 2 


- = i (tan e r + tan e A ) 


X "B 


yields 


X* °X* 


Since this flow is nonuniform a characteristic calculation similar to an 
interior point calculation yields the flow properties at C^. Note that 
point E on the C + characteristic EC^ is interpolated between B and A^. 

The jump relations (Equations 38 - 43) are solved using the determined up- 
stream conditions based on the assumed angle This yields all properties 

at C£. Using the deflection angle calculated from the jump conditions, 

a C + characteristic calculation performed along (F-Cp) yields an alternate 
value of the pressure at p^ . The pressures are compared and if the dif- 
ference exceeds a specified tolerance, a new value of is assumed and 
the- process repeated until convergence is obtained. After convergence with 
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1 1 

a = 1, a = 0, the calculation is repeated using a = 6 * 

C. Under-Expansion Interaction - The program developed has the 
capability for equilibrium flows of computing the under-expansion interaction 
produced by pressure mismatch between the nozzle and a surrounding airstream. 
This situation is depicted in Figure (8a). Under-expansion conditions occur 



as a result of either P. > or e. > e or some combination of both condi- 

j e j e 

tions. Generally P- > > P defines an under-expanded flow. It is assumed 

j e 

that during the under-expansion interaction, the species remain chemically 
in equilibrium. The expansion is isentropic and the local interaction is 
two dimensional and inviscid in the limit of vanishing radial distance with 
respect to the cowl edge. 

The basic equations describing the Prandtl-Meyer expansion process are 

p/p r - constant ( 44 ) 
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h + 1 ? = constant (45) 

^ + \ d(V 2 ) = 0 (46) 

i d ln(p) ± d 6 = 0 (47) 

r = r (p, h,$) where $ = constant (14) 


For a small incremental step Ap, Equations (46), (47) 14) and (50) can be 
written 


- V 1 = - <£) 2 -£V 

(48) 

jT In (p£/P] > ) - ( ® 2 _ ® i ) = 0 

(49) 

P2^ p 2 ~ P l^ p l 

(50) 


where r is held constant in the integration step, yielding values for V^i 
pg and 0g, Then Equation (45) yields h 9 . In this manner, the Prandtl- 
Meyer equations may be integrated, for a mixture in equilibrium. 

Since the flow deflection and pressure downstream of the shock wave and 
Prandtl -Meyer are unknown an iteration process is required. A typical in- 
teraction calculation proceeds as follows. A shock wave angle is assumed 
for which flow properties (p, h, e etc.) are computed downstream of the shock 
wave. Equations (44) through (47) are solved using small increments of Ap. 
The pressure behind the shock is the final pressure and the jet pressure is 
the initial pressure. If the turning angle for the expansion does not agree 
with the flow angle behind the shock to within a specified tolerance, a new 
shock wave angle is assumed and the process repeated until convergence is 
obtained. 
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After this solution is obtained, the program continues the normal calcula- 
tion procedure until the last expansion ray is completed. The program 
terminates the calculation along the last ray. It is anticipated that the 
flow along the dividing streamline will not in general affect the pressure 
distribution along the undersurface, for the vehicles to be considered. 

D. Over-expansion Phenomena - The nozzle over-expansion at the 
cowl is computed in a similar fashion to the under-expansion phenomena, 
except that a shock wave is required in the nozzle flow and an expansion in 
the external flow as depicted in Figure (8b). For purposes of simplicity 



FIGURE 8b. OVER-EXPANSION INTERACTION 

it is assumed that the external flow is initially uniform with constant ratio 
of specific heats. Further, it is assumed that pressure-flow deflection re- 
lationship on the external side of the dividing streamline is described by 
simple Prandtl-Meyer relations. These assumptions do not inhibit the programs 
generality, but are intended only to simplify the computation. They are 
readily removed and more general but complex methods can be used if the need 
arises . 
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Since the program stores data on C_ characteristics, shock waves of the C + 
family are conveniently traced Thus, in order not to disrupt the program 
logic the problem is inverted so as to trace the over-expansion shock as a 
C, wave. This procedure is performed automatically as part of the program 

T 

logic. 


E. Contact Surface - A contact surface is a stream surface of 
the flow, therefore, the pressure and flow deflection must be continuous 
across the discontinuity. Figure (9) illustrates a contact surface calcula- 



FIGURE 9. CONTACT CALCULATION 

tion for supersonic flow. In the present program the characteristic on the 
external flow side is replaced by a Prandtl-Meyer pressure-flow deflection 
relation and the external flow is assumed uniform. The solution requires a 
iterative procedure similar to a wall boundary calculation except that the 
shape of the boundary (i .e. , contact surface) is not known a priori. 
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In Figure (9) CD is the contact and point A lies on a C_ characteristic 
CA, on which all properties have previously been calculated. A guess is 
made for 0p and for a = 1 and M 0 a simultaneous solution of the equations 


and 


^ = I {tan S D + tan 0 C ) = M 3 


y D~ y A 

= M. 

X D‘ X A 1 nozzle 


( 51 ) 


(18a) 


yields XD and YD. Using the guessed value of e D from the characteristic re- 
lation along AC and the Prandtl -Meyer relation in the external stream two 
values of pressure are obtained at D. If these do not agree to within a 
specified tolerance a new guess for is made and the process repeated ■ 

until convergence is achieved. Using the streamline relations along CD the 
remaining properties (q, h, T etc.) are computed on each side of the dis- 
continuity and the process repeated for a = p 3 = 


F. Shock Reflection At Wall - The incident and reflected strength 
of a shock wave at a wall boundary is determined by the condition that down- 
stream of the reflected wave (3) the flow deflection at the wall must equal 
the wall slope; Figure (10) depicts this interaction. 



SHOCK CCO 


FIGURE 10. SHOCK REFLECTION AT WALL 
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The solution requires an iterative procedure since the reflected wave 
strength is a function of the data downstream of the incident wave (2) and 
subject to the above constraint. A shock angle sigma (c^) is assumed for 
the incident wave and the location of the shock wall intersection (x D » y D ) 
is computed using 


V y A _ 1 


= (tan + tan c^) = Mg 


X D~ X A "1 


and 


y D = A j ^ x D' x i^ + B j ^ X D' X 1^ + C j 


which yields (Xp, y D , 6p). 


Since the flow upstream of the shock wave is nonuniform, a characteristic 
solution similar to a wall calculation is required to determine the flow 
properties at D^. Note that point A is interpolated on the C_ character- 
istic CB^. Then the Hugoniot relations Equations (38) through (43) yield 
the flow properties downstream of the incident wave (2). The pressure from 
this calculation is then compared with a characteristic calculation along 
ED 2 on the downstream side. Point E is interpolated on C_ characteristic 
B^F , If the pressures do not agree to within a specified tolerance, a new 
shock angle is assumed and the process repeated until convergence. 


These properties are then used as upstream conditions for the reflected wave. 
Assuming a reflected wave angle sigma (og) the Hugoniot relations yield the 
flow properties downstream. If the flow angle does not agree with the wall 
angle to within a specified tolerance a new reflected wave angle is assumed 
and the process repeated until convergence. After convergence the entire 
calculation is repeated with a « \ and B = j- 

Since the program logic stores data only on C_ characteristics, shock waves 
of the C + family are conveniently traced. However, the reflected wave is a 
C wave as depicted in Figures (3) and (10). In order not to disrupt the 
program logic, after computing the strength of the reflected wave the flow 
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field is automatically inverted making the reflected wave a C + shock wave. 
Special provisions are required, however, for storing new initial data along 
a line which enables the propagation strength of the reflected wave to be 
computed up to the next boundary. The program logic is such that these re- 
quirements are performed automatically. The overall grid for this calcula- 
tion is shown in Figure (3) of Section IIC. 
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IV. SAMPLE CALCULATIONS 


Figure (11) depicts an under-expanded nozzle calculation using NASA sup- 
plied geometry and initial conditions. The upper curves represent a com- 
parison of vehicle undersurface pressures between two dimensional NASA 
calculations and the subject ATL program. Excellent agreement is seen 
to exist over the length of the undersurface. 

The lower curve in Figure (11) represents a nozzle with the same geometry 
and initial conditions except that a lateral nozzle area variation has 
been provided for between the throat surface and the cowl. The origin of 
the source is at x/h^ - 7. Downstream of the cowl the flow is assumed to 
be two dimensional. The lateral area variation is seen to produce sig- 
nificant changes in the pressure distribution on the undersurface and in 
the location and strength of waves. Thus, this approximation to the 
lateral area variation can be a powerful tool in the design of scramjet 
nozzle exhaust flow fields. 

Figures (12) and (13) demonstrates the programs capability to calculate 
over-expanded nozzle flows. Figure (12) is a trace of the vehicle geometry, 
shock shape and contact shape upto x/h^ = 27. Figure (13) indicates the 
pressure distributions on the cowl - contact surface and the vehicle under- 
surface for this case. The under-expansion shock jump is clearly visible 
at the cowl trailing edge (x/h t =6). 
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A. Case la - Two Dimensional 


Initial Conditions 



Free Stream Conditions 

Temperature 

(T/TJ 

- 

10.1 

Altitude 

101,800 ft. 

Pressure 

(P/Pj 

- 

36.65 

Mach Number M 
00 

10. 

Velocity 

(q/uj 

- 

0.929 

Pressure p 

r CO 

23.09 lb/ft 

Mach Number 

M 

- 

2.91 

Temperature - 

418. 8°R 

Flow Angle 

e 

- 

0. 



Equivalence Ratio <i> 

- 

1.0 



Source origin 

- none* 

two dimensional 



External Flow 






Temperature 

(T/TJ 

- 

1.0 



Pressure 

(P/P.) 

- 

1.0 



Velocity 

(q/O 

- 

1.0 



Mach Number 

M 

- 

10. 



Flow Angle 

e 

- 

0. 



Equivalence Ratio $ 

- 

0. 



Geometry - Y=AX 2 +BX+C 






Vehicl 

e 


Cowl 


"vf 

. 

1 

O 

X 

.4 - 8.0 


8.0 - End. 

0 - .4 .4 - : 

3.0 

A **.5565 

0 


.01019 

.1314 0 


B 0 

-.4452 


-.6082 

0 .1051 


C 0 

.08905 


.7410 

1.0 .9790 



B. Case lb - Source Flow 

All data is the same as Case (la) except initial source origin is located 
at x = -7.0. 



Wall Pressure, psf 
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Case II - Two Dimensional 


Initial Conditions 

Free Stream Conditions 

Temperature 

(T/TJ - 10.77 

Altitude 

- 60,000 ft. 

Pressure 

(P/P J - 20.55 

Mach Number 

- 4.0 

Velocity 

(q/uj - -904 

Pressure p 
1 00 

- 151 Ibs/ff 

Mach Number 

M - 1.09 

Temperature T m 

- 390°R 

Flow Angle 

S - 0. 



Equivalence Ratio $ - 0. 



Source origin 

- none, two dimensional 



External Flow 




Temperature 

(T/TJ - 1.327 



Pressure 

(P/P J - 2.493 



Velocity 

(q/uj - -95 



Mach Number 

M - 6.67 



Flow Angle 

0 - 0. 



Equivalence Ratio $ - 0. 



Geometry - Y=AX^+BX+C 




Vehicle 

Cowl 


X 0 - ,4 

.4 - 8.0 8.0 - End 

0 - . 4 

.4 - 6.0 

A -.5565 

0 .01019 

.3117 

0 

B 0 

-.4452 -.6082 

0 

.2493 

C 0 

.08905 .7410 

1.0 

.9501 
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V. CONCLUSIONS 

A unique characteristic procedure has been developed which computes realistic 
scramjet nozzle exhaust flow fields. The use of an axially expanding co- 
ordinate system (line source) allows lateral nozzle area variations to be 
accounted for in a quasi two dimensional fashion. Additional geometric 
flexibility is incorporated in the numerical procedure through the use of 
additional coordinate systems i.e., axisymmetric, and Cartesian. 

The numerical procedure uses a free running characteristic grid, but stores 
data on the C_ characteristics. The technique is second order accurate in 
the characteristic sense. 

The wide variety of boundary conditions incorporated into the program permits 
the calculation of wall boundaries, shock boundaries, contact boundaries, 
shock-wall intersections and over or under-expansion interactions. 

The use of the unique hydrogen-air equilibrium curve fits developed in Refer- 
ence (5), as well as frozen hydrogen-air chemistry permits bounding the com- 
plex chemical phenomena which occur in exhaust nozzle flows. 

Thus, it is felt that the current program will give the designer of scramjet 
nozzle exhaust flow fields a flexibility not available with previous methods. 
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APPENDIX I . • 

CHARACTERISTIC DERIVATION , 

FROZEN FLOW ' 

The continuity equation (Equation (1)) may be expanding yield- 
ing: 

+ J M S i n 0 T* J 2 ^ cos 9 = 0 (1.1) 


» U * - H * « H 


From the s-momentum equation 


0 . Il£ 

p 3S q • 3 s 


and using both the Equation of State (6) and the energy Equa- 
tion (4), the term may be replaced by - 


-3p 

3 S 


W Y, 

TT 


(■Vco-l) 

~ 



Making these substitutions in Equation (1.1), and multiplying 
through by q, we obtain 


W Y 


„ 


(y:.-D 


se 


M__ - 1 ) l£ + p q ~ £_E_ 

' 1 3S - 8n 


(1.2) 


1 2 
+ ~ q 2 sin 9 + J 2 cos e = 0 


By algebraic manipulation, using Equations (11-) and (12), the 
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A ElLENDI X I (Continued! 


term multiplying |-jr in (1.2) may be reduced to (M 2 -!) and 
hence Equation (1.2) becomes 


(M 2 -l) 3 p . 

, 1_ 4. 

P q^ 9s 


a e 
an 


+ 


*s i n 



0 


Using Equations (11) and (14) pq 2 = ypM 2 . . Making this sub- 
stitution in both Equations (1.3) and the normal momentum 
equation, we obtain 


(1*3) 



(M 2 -!) 

YpM 2 



ae 

an 



0 - 



0 




ii - 
as 



The total derivatives of p and e may be expressed by: 


i£, 

ds 

+ iR. 

dn 

= dp 

as 


3n 


39 

d$ 

+ i£ 

dn 

= de 

as 


an 




(1.4) 


(1.5) 


( 1 . 6 ) 

(1.7) 


Written in matrix form, the above system (Equations (1.4) through 
(1.7)) becomes 
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APPENDIX I (Continued) 

\ 4 

* 

j( M 2 -l )/ypM 2 0 

0 1/yM 2 P 

ds dn 




(i:s) 


The characteristic directions of this system of equations are 
obtained by setting the determinant of the coefficient matrix 
equal to zero , 


(M 2 -1)/yM 2 p 0 01 

0 1/yM 2 p 1 0. 

t 

ds dn 0 0 


= 0 ■ - (1.9) 


0 


0 


ds dn 
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dn^ 

obtaining “ 


(M 2 -l) 


whose slope is given by 


Hence the characteristics are lines 


dn s 
ds 




tann 


or expressed in Cartesian coordinates 


a tan ( e±y ) 

The compatibility relation along the characteristics is ob- 
tained by replacing any column of the coefficient matrix with 
the vector on the right-hand side of Equation (1.8) and setting 
the determinant of this matrix equal to zero. 


0 . s i n e 

-<-* 7 — 0 0 1 

°2 

+ cose) 

A 


0 

dp 


dn 0 0 


= 0 


I de 0 ds dn | 

Expanding the determinant and using Equation (1.10) we obtain 
the compatibility relation (Equation 8). 


J,sine J 0 cose 

S l. n» c o s . y d ln p ± d6 + ( 1 - + ^ — ) 


i — SillH — ^ d x = 0 
(cose±u ) 


( 1 . 10 ) 


( 1 . 11 ) 


( 1 . 12 ) 


(1.13) 
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Since we may write total derivatives for the changes in entropy, 
stagnation enthalpy and species mass fraction along streamlines, 
the streamlines act as characteristics in the flow field. Hence, 
along a streamline, whose slope is given by 


P- - tane ’ (1.14) 

dx 

the following equations hold: 

dS = 0 (1.15) 

dH = 0 (1.16) 

and da.j = 0 ( 1 . 17a ) 


For a flow in chemical equilibrium the derivation is identical 
with the frozen ratio of specific heats y replaced by r and the 
Mach number defined in terms of the local equilibrium sound 
speed. Equation (1. 17a) is replaced by 


d$> = 0. 


(1 . 17b) 
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APPENDIX I la 


THERMODYNAMIC COEFFICIENTS ' 

F.OR FROZEN SPECIES 

The following thermodynamic properties have been tabulated as 
polynomials in temperature (in degrees Kelvin) in the form be- 
low: 

C* 

= a 1 + a 2 T* + a 3 T* 2 + a 4 T* 3 + a 5 T* 4 (Ila.-l) 

R 0 


RnT' 


= a 


. tl T* + f-3 t *2 + i4 t *3 + ii t *4 + *6 (Ila. 2) 

1 2 3 4 5 - j* 


a 1 In T* + a 2 T* 


+ ll T* 2 + — T* 3 


(Ila. 3) 


+ £5 T* 4 + a 7 

A / . 


The coeffients (a - ! “ a^) have been tabulated for the temperature 
intervals 300° to 1 000°k and 1 000° to 5000°k in Reference (4). 
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• APPENDIX 1 1 b 
EQUILIBRIUM HYDROGEN -AIR 
CURVE FITS FOR f‘, ' h and p - 


The variation of T (the equilibrium value of y) as a func- 
tion of temperature (T), pressure (P) and equivalence ratio 
($) is presented graphically in Figures (A I ) , (A2 ) and (A3) 
from values tabulated in Reference (4). In Figure (A1 ) it 
can be seen that r is a strong function of T over the tem- 
perature range of interest, while the effect of varying com- 
position is small by comparison. Moreover, Figure (A2.) in- 
dicates that r is moderately sensitive to pressure and the 
degree of sensitivity increases substantially as the tempera- 
ture level increases and dissociation effects become impor- 
tant. 


As a result of these observations, temperature is the pri- 
mary. independent variable, while pressure is the secondary 
independent variable and composition acts as a perturbation 
variable. Thus, we can fit the function r (T,P,$) with a 
polynomial in T and add on a temperature dependent correc- 
tion term for the effect of pressure and a temperature in- 
dependent correction term for the effect of $ . 

An examination of Figure ( A 1 ) suggest that the function T (T) 
can best be curve fit by breaking up the temperature range 
into three intervals such that the function can be represent- 
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ed by a parabola in each range. Choosing p - IQ 5 pascal and 
0=1 as our base, we therefore find three functions 

„ ‘ i 

r 1 (T,10 5 ,l) = - 1.833xl0“ 7 T 2 + 7.5xl0" 5 T + 1.367 

“■ 4- 


r 2 (T,10 5 ,l) = 2.0xl0“ 8 T 2 - 1.38X10- 4 T + 1.423 
r 3 (T,10 5 ,l) = 7 . 27x10 “ 8 T 2 - 4.57xl0“ 4 T + 1.85 


and define the basic temperature function as 


r(T,io 5 ,D =< 


Tid.l 0 5 ,1) 
r 2 (T,10 5 ,l) > for 
r 3 (T ,io 5 ,1) 


' T«500°K 

V 500<T<2000°K > 
T>2000°K , 


Figure (A3 } indicates that 3£ is constant in the two ranges 
■ - 3 $> 

$<1 and.$>l, but is a function of T. Fitting the function 

: ar in each of the ranges of <f> we obtain 

aT 


*n i ( T )" 


r "t 

$£l 

for 


_ n 2 ( T )_ 


0>1 
L — J 


where 

n (T) = 4xlO" 9 T 2 - 2xl0" 5 T - 0.019 
„ (T) = 3 . 39x 10" 2 T°- 5 - 3.91xl0' 4 T -0.681 


(Ilb.l) 
(IIb.2) 
(lib. 3) 

(IIb.4) 


(iiti. 5) 

(lib. 6) 
(IIb.7) 
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This now defines r as a function of both temperature and $ 

by means of the equation ' 

* ' • ' 

r(T , 10 5 ,4>) = r(T,10 5 ,l) + Oil) §| {ll'b.8) 

Finally, the effect of pressure must be included. From 
Figure "(18) we observe that r may be approximated as 

■ rCT.P.b) = r(T,10 5 ,*) •+ m tlo^(p) -5] {lib. 9) 

where m is a function of T. Deriving m , we find 

|t<iooo°k 

for . . _ r (Mb. 10) 
T>1000°K 

Summarizing, the final function obtained is 

- - • r(T,P,«) = r(T,10 5 ,l) + m(l|-P -5) + || («-l) 

f 

, i j- 

wherethe functions r(T,10 5 ,l), 3T_ and m are given by Equa- 
tions '(4), (5) and (10) respectively. 

1 J * 

The curve fit for enthalpy is derived in a similar way. Figures 

(A-4) and (A-5) present the variation of h with temperature, 
pressure and equivalence ratio. As was the case for r» the 


0 

-2. 15x1 0"M + 0 .9 1x10" 4 T- 0.0695 




FIGURE A4 . ENTHALPY AS A FUNCTION OF TEMPERATURE (p=10 5 pa.) 
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PRESSURE (pa) 

FIGURE A5 , ENTHALPY AS A FUNCTION OF PRESSURE. 
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function h(T, p) is fit by a quadratic function of T, the 

coefficients of which are functions of 4> and an additive term 

* / 

for the effects of pressure. The resulting curve fit is sum- 
marized below, • ' 


h(T,<j>,p) 


where h(T,<i>,p) = h ( T A, IQ 5 ) 1 + 


h(T,<f>,10 5 ) 


T- 2000°K 

< 

► for < 


h(T.<f,p) . 


n 

T>2000°K 


,( l + ft ) ( T- 2000 ) 
2000 


■ 125( 2T3 -5) 2 -. 275(^|-5)| 


The basic function hfT^lO 5 ) is defined as 

h (T ,4> , 10 5 ) = 10 6 (a,T 2 + b i T + c 1 ) 

with the coefficients a ly bj and c 1 defined below: 
For ' T * 2000°K and <p < 1 


a l = 

1CT 7 (-. 1042<f> 2 + 

.82424) +.987) 

b r - 

10‘ 3 ( .01167<J> 2 + 

. 1503$ + .938) 

c l- = 

. 0284<f. 2 + . .6731 

<(> + .4293 


(lib. 11) 


(lib. 12) 


(lib. 13) 


(lib. 14) 


(lib. 15) 
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For T < 2000°K and <j> > 1 

a j = ' 10' 7 { 1 . 7874> 2 - 5. 48$ + 5.4) 
bi « 10" 3 (-.1867$ 2 + 1.11$ + .176) 

c. = -.0933$ 2 + 3.975$ - 2.808 

For T > 2000°K and $ i 1 

a, = 10' 6 (1.792$ 2 + .3983$ + .310) 
b x = 10’ 3 ( -9 . 05$ 2 - .07917$ + .245) 
c, = 10.86$ 2 - .1183$ + .970 

For ’ T > 2000°K and $ >1 - 

aj = 10' 6 ( 4 . 81$ 2 - 13.9$ +11.59) 

" • . bj = 10“ 3 (-23.08(j) 2 + 66.8 2<J> - 52.61) 

Cl = 27 . 054> 2 - 73.73^ + 58.39 


(lib. 16) 


(rib'. 17) 


(lib. 18) 


When the inverse function T(h # 4>,p) is required, it is obtained 
by an iterative solution of Equations (12) through (18) 


The density is found by obtaining a curve fit for the mixture 
molecular weight and using the equation of state 
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P 




(lib. 19) 


RT 




Where R is the universal gas constant and m is the 

molecular 

weight. 

• '■ 

.» ' 



The behavior of m with 

T,p and <(> is illustrated in 

Figures (A6) 

and ( A7 ) . We see that 

for temperatures less 

than 

2000°K, 

m i s 

essentially independent of temperature. The 

di scon ti nu i tv i n 

slope of m(4>) shown in Figure (A7) requires that the equivalence 
ratio range be split in two. Thus, 

for T s 2000°K 





• m(4>) = < 

| 1.53cj> 2 -5.895<i> + 28.965 


U t i 


1.60cJ> 2 -10.6<j) + 33.6 

► for 

♦ > 1 

* (lib. 20) 

For thejiigher temperature range, it is convenient 
the form 

to empl oy 

- ■ / 

/* m = 

(S( p . <f> ,T ) 



'(lib. 21) 

where 

t 





_ n 0 (4>) 




6 » d 2 (p.*) ( T “ 2000 ) ‘ 

1 K 1000' 



(lib. 22) 







«b 


FIGURE A7 . MOLECULAR WEIGHT AS A FUNCTION OF EQUIVALENCE RATIO FOR T<2000°K. 
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and 

• d2 = a 2 c^rf) 1 ' 5 + b 2 (|t§) + c 2 


: : a 2 

b 2 

For 0 £ 4 csl c 

■ c 2 

n 2 


® 2 

b 2 

and for 1 < 4<2 

c 2 

- n 2 


- 2 . 3 <J> 2 + 4 . 01 d> + 1 . 736 
8 . 6 1 d> 2 . 15 . 4 24 > - 6.66 
- 16 . 88 c )) 2 + 33 . 21 c)) + 14.58 
. 43754 2 + .06254 + 2.08 

- . 3224 2 + 2.3634 + 1.905 
2 . 764 2 -. 7.564 - 8.68 
3 . 64 ^ + 7.364 + 27.15 
. 474 2 + 1 .8254 + .350 


h 

i 


14 


(lib. 23 } 


(lib. 24 ) 


(lib. 25 ) 
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APPENDIX III 

THRUST, LIFT AND PITCHING MOMENT 


In general, the thrust, lift and pitching moment can be defined from the follow 

ing pressure integral taken over all the nozzle surfaces up. to the final X sta- 
tion (XFINAL) , i . e . : 




V dA n 


(III. la) 


L y = 


(P-Pj y dA n 


(III. lb) 




(P‘Pj i y ‘XdA n + 


(p-pj i x -ydA n 


(III.lc) 


i 

where the coordinate system and vehicle configuration are depicted in Figure 
(III-l). 


However, since the lateral geometry of the nozzle is treated approximately, it 
is not possible to determine the lateral ..contributions to the aboVe’integrals by 
direct pressure integration. The use of the integral conservation theorems 
provides, however, an alternate means of defining the thrust, lift and pitching 
moment. 


For a fixed control volume 






(q x r) ( q)-ndA 


(III. 2a) 


ifx? 


(III. 2b) 
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where the integrals extend over the throat area and a suitably defined noz- 
zle exit area. Care must be exercised in using Equations (III. 2} since small 
errors in mass flow can produce large errors in net thrust, lift and pitching 
moment. 


By straightforward algebraic manipulation and use of the Equation of State and 
the definition of the local sound speed Equations (III. 2) may be reduced to 
the following: 


T x = 


2 

(ypM sin(0 c -e) cose + (p-p ) sine e ) 

^ o - 


sine. 


dydz 


^exit 


(ill. 3a) 


2 

(vpM sin(e-e) cose + (p-p ) sine ) 

5 _________ 00 _5 

sTnT"^~ 


dydz 


«A 


throat 


L y = 


2 

(ypM sin(e c -e) sine + (p-p ) cose ) 

oo 5 


sine. 


dydz 


exit 


(III. 3b) 


2 

(ypM sin(e c -e) sine + (p-p ) cose_) 

oo 5 


A 


sine. 


dydz 


throat 
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where e g is the local inclination of the throat or exit area and z is the 
lateral extent of the nozzle. The algebraic differences between Equations (111.2) 
and (I II. 3) represent the integrated force and moment contributions of the side- 
walls/and or fences. If no fences are present the momentum balance is carried 
out at the cowl exit station but the pressure integrations are still computed 
over the full vehicle and cowl surfaces. 




VEHICLE 

UNDERSURFACE 




FIGURE III-l. THRUST, LIFT, MOMENT 
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